{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TLncjuFUlkce"
   },
   "source": [
    "**Install deepxde**  \n",
    "Tensorflow and all other dependencies are already installed in Colab terminals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5635,
     "status": "ok",
     "timestamp": 1640283886921,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "dibZzRjplLH-",
    "outputId": "226fcb6f-10e7-4637-e96e-f20e0a6e5b5d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting deepxde\n",
      "  Downloading DeepXDE-0.14.0-py3-none-any.whl (111 kB)\n",
      "\u001b[?25l\r\n",
      "\u001b[K     |███                             | 10 kB 22.0 MB/s eta 0:00:01\r\n",
      "\u001b[K     |█████▉                          | 20 kB 21.5 MB/s eta 0:00:01\r\n",
      "\u001b[K     |████████▉                       | 30 kB 10.9 MB/s eta 0:00:01\r\n",
      "\u001b[K     |███████████▊                    | 40 kB 8.5 MB/s eta 0:00:01\r\n",
      "\u001b[K     |██████████████▊                 | 51 kB 5.5 MB/s eta 0:00:01\r\n",
      "\u001b[K     |█████████████████▋              | 61 kB 5.6 MB/s eta 0:00:01\r\n",
      "\u001b[K     |████████████████████▋           | 71 kB 5.5 MB/s eta 0:00:01\r\n",
      "\u001b[K     |███████████████████████▌        | 81 kB 6.2 MB/s eta 0:00:01\r\n",
      "\u001b[K     |██████████████████████████▍     | 92 kB 4.9 MB/s eta 0:00:01\r\n",
      "\u001b[K     |█████████████████████████████▍  | 102 kB 5.3 MB/s eta 0:00:01\r\n",
      "\u001b[K     |████████████████████████████████| 111 kB 5.3 MB/s \n",
      "\u001b[?25hRequirement already satisfied: scikit-learn in /usr/local/lib/python3.7/dist-packages (from deepxde) (1.0.1)\n",
      "Requirement already satisfied: matplotlib in /usr/local/lib/python3.7/dist-packages (from deepxde) (3.2.2)\n",
      "Collecting scikit-optimize\n",
      "  Downloading scikit_optimize-0.9.0-py2.py3-none-any.whl (100 kB)\n",
      "\u001b[K     |████████████████████████████████| 100 kB 8.9 MB/s \n",
      "\u001b[?25hRequirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from deepxde) (1.19.5)\n",
      "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from deepxde) (1.4.1)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->deepxde) (3.0.6)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->deepxde) (2.8.2)\n",
      "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib->deepxde) (0.11.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->deepxde) (1.3.2)\n",
      "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.1->matplotlib->deepxde) (1.15.0)\n",
      "Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.7/dist-packages (from scikit-learn->deepxde) (1.1.0)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.7/dist-packages (from scikit-learn->deepxde) (3.0.0)\n",
      "Collecting pyaml>=16.9\n",
      "  Downloading pyaml-21.10.1-py2.py3-none-any.whl (24 kB)\n",
      "Requirement already satisfied: PyYAML in /usr/local/lib/python3.7/dist-packages (from pyaml>=16.9->scikit-optimize->deepxde) (3.13)\n",
      "Installing collected packages: pyaml, scikit-optimize, deepxde\n",
      "Successfully installed deepxde-0.14.0 pyaml-21.10.1 scikit-optimize-0.9.0\n"
     ]
    }
   ],
   "source": [
    "!pip install deepxde"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H9h3CaW8kj1y"
   },
   "source": [
    "**Problem setup**  \n",
    "  \n",
    "We are going to solve the non-linear Schrödinger equation given by  \n",
    "$i h_t + \\frac{1}{2} h_{xx} + |h|^2h = 0$  \n",
    "  \n",
    "with periodic boundary conditions as  \n",
    "$x \\in [-5,5], \\quad t \\in [0, \\pi/2]$  \n",
    "$h(t, -5) = h(t,5)$  \n",
    "$h_x(t, -5) = h_x(t,5)$  \n",
    "  \n",
    "and initial condition equal to  \n",
    "$h(0,x) = 2 sech(x)$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f7Fo1dqAmMt-"
   },
   "source": [
    "Deepxde only uses real numbers, so we need to explicitly split the real and imaginary parts of the complex PDE.  \n",
    "  \n",
    "In place of the single residual  \n",
    "$f = ih_t + \\frac{1}{2} h_{xx} +|h|^2 h$  \n",
    "  \n",
    "we get the two (real valued) residuals  \n",
    "$f_{\\mathcal{R}} = u_t + \\frac{1}{2} v_{xx} + (u^2 + v^2)v$  \n",
    "$f_{\\mathcal{I}} = v_t - \\frac{1}{2} u_{xx} - (u^2 + v^2)u$  \n",
    "  \n",
    "where u(x,t) and v(x,t) denote respectively the real and the imaginary part of h.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 6652,
     "status": "ok",
     "timestamp": 1640283904989,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "YjOFqj77nl_R",
    "outputId": "aa0997be-fdbe-4a21-8def-580cf5e20c6c"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Deepxde backend not selected or invalid. Assuming tensorflow.compat.v1 for now.\n",
      "Using backend: tensorflow.compat.v1\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setting the default backend to \"tensorflow.compat.v1\". You can change it in the ~/.deepxde/config.json file or export the DDEBACKEND environment variable. Valid options are: tensorflow.compat.v1, tensorflow, pytorch (all lowercase)\n",
      "WARNING:tensorflow:From /usr/local/lib/python3.7/dist-packages/tensorflow/python/compat/v2_compat.py:111: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "non-resource variables are not supported in the long term\n",
      "WARNING:tensorflow:From /usr/local/lib/python3.7/dist-packages/deepxde/nn/initializers.py:116: The name tf.keras.initializers.he_normal is deprecated. Please use tf.compat.v1.keras.initializers.he_normal instead.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import deepxde as dde\n",
    "\n",
    "# For plotting\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import griddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "executionInfo": {
     "elapsed": 337,
     "status": "ok",
     "timestamp": 1640283910535,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "Oj-6F7RuobUn"
   },
   "outputs": [],
   "source": [
    "x_lower = -5\n",
    "x_upper = 5\n",
    "t_lower = 0\n",
    "t_upper = np.pi / 2\n",
    "\n",
    "# Creation of the 2D domain (for plotting and input)\n",
    "x = np.linspace(x_lower, x_upper, 256)\n",
    "t = np.linspace(t_lower, t_upper, 201)\n",
    "X, T = np.meshgrid(x, t)\n",
    "\n",
    "# The whole domain flattened\n",
    "X_star = np.hstack((X.flatten()[:, None], T.flatten()[:, None]))\n",
    "\n",
    "# Space and time domains/geometry (for the deepxde model)\n",
    "space_domain = dde.geometry.Interval(x_lower, x_upper)\n",
    "time_domain = dde.geometry.TimeDomain(t_lower, t_upper)\n",
    "geomtime = dde.geometry.GeometryXTime(space_domain, time_domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "executionInfo": {
     "elapsed": 283,
     "status": "ok",
     "timestamp": 1640283914251,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "dbx1YulhoORs"
   },
   "outputs": [],
   "source": [
    "# The \"physics-informed\" part of the loss\n",
    "\n",
    "\n",
    "def pde(x, y):\n",
    "    \"\"\"\n",
    "    INPUTS:\n",
    "        x: x[:,0] is x-coordinate\n",
    "           x[:,1] is t-coordinate\n",
    "        y: Network output, in this case:\n",
    "            y[:,0] is u(x,t) the real part\n",
    "            y[:,1] is v(x,t) the imaginary part\n",
    "    OUTPUT:\n",
    "        The pde in standard form i.e. something that must be zero\n",
    "    \"\"\"\n",
    "\n",
    "    u = y[:, 0:1]\n",
    "    v = y[:, 1:2]\n",
    "\n",
    "    # In 'jacobian', i is the output component and j is the input component\n",
    "    u_t = dde.grad.jacobian(y, x, i=0, j=1)\n",
    "    v_t = dde.grad.jacobian(y, x, i=1, j=1)\n",
    "\n",
    "    u_x = dde.grad.jacobian(y, x, i=0, j=0)\n",
    "    v_x = dde.grad.jacobian(y, x, i=1, j=0)\n",
    "\n",
    "    # In 'hessian', i and j are both input components. (The Hessian could be in principle something like d^2y/dxdt, d^2y/d^2x etc)\n",
    "    # The output component is selected by \"component\"\n",
    "    u_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)\n",
    "    v_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)\n",
    "\n",
    "    f_u = u_t + 0.5 * v_xx + (u ** 2 + v ** 2) * v\n",
    "    f_v = v_t - 0.5 * u_xx - (u ** 2 + v ** 2) * u\n",
    "\n",
    "    return [f_u, f_v]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "executionInfo": {
     "elapsed": 323,
     "status": "ok",
     "timestamp": 1640283979042,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "OisdDqf0oSLx"
   },
   "outputs": [],
   "source": [
    "# Boundary and Initial conditions\n",
    "\n",
    "# Periodic Boundary conditions\n",
    "bc_u_0 = dde.icbc.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0\n",
    ")\n",
    "bc_u_1 = dde.icbc.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0\n",
    ")\n",
    "bc_v_0 = dde.icbc.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=1\n",
    ")\n",
    "bc_v_1 = dde.icbc.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=1\n",
    ")\n",
    "\n",
    "# Initial conditions\n",
    "def init_cond_u(x):\n",
    "    \"2 sech(x)\"\n",
    "    return 2 / np.cosh(x[:, 0:1])\n",
    "\n",
    "\n",
    "def init_cond_v(x):\n",
    "    return 0\n",
    "\n",
    "\n",
    "ic_u = dde.icbc.IC(geomtime, init_cond_u, lambda _, on_initial: on_initial, component=0)\n",
    "ic_v = dde.icbc.IC(geomtime, init_cond_v, lambda _, on_initial: on_initial, component=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "executionInfo": {
     "elapsed": 259,
     "status": "ok",
     "timestamp": 1640283983373,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "u-f8sfj2oWgy"
   },
   "outputs": [],
   "source": [
    "data = dde.data.TimePDE(\n",
    "    geomtime,\n",
    "    pde,\n",
    "    [bc_u_0, bc_u_1, bc_v_0, bc_v_1, ic_u, ic_v],\n",
    "    num_domain=10000,\n",
    "    num_boundary=20,\n",
    "    num_initial=200,\n",
    "    train_distribution=\"pseudo\",\n",
    ")\n",
    "\n",
    "# Network architecture\n",
    "net = dde.nn.FNN([2] + [100] * 4 + [2], \"tanh\", \"Glorot normal\")\n",
    "\n",
    "model = dde.Model(data, net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hu5Tyq32DezT"
   },
   "source": [
    "Adam optimization.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 439512,
     "status": "ok",
     "timestamp": 1640284427354,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "BR_s0D-_oZYz",
    "outputId": "9071ee58-c123-427b-fe32-cc8fa7c721f4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model...\n",
      "Building feed-forward neural network...\n",
      "'build' took 0.103608 s\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/dist-packages/deepxde/nn/tensorflow_compat_v1/fnn.py:110: UserWarning: `tf.layers.dense` is deprecated and will be removed in a future version. Please use `tf.keras.layers.Dense` instead.\n",
      "  kernel_constraint=self.kernel_constraint,\n",
      "/usr/local/lib/python3.7/dist-packages/keras/legacy_tf_layers/core.py:255: UserWarning: `layer.apply` is deprecated and will be removed in a future version. Please use `layer.__call__` method instead.\n",
      "  return layer.apply(inputs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'compile' took 4.128636 s\n",
      "\n",
      "Initializing variables...\n",
      "Training model...\n",
      "\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\n",
      "0         [1.87e-02, 1.38e-03, 5.78e-03, 4.22e-05, 9.99e-01, 1.20e-04, 7.86e-01, 1.25e-01]    [1.87e-02, 1.38e-03, 5.78e-03, 4.22e-05, 9.99e-01, 1.20e-04, 7.86e-01, 1.25e-01]    []  \n",
      "1000      [1.03e-02, 1.17e-02, 2.13e-04, 7.70e-05, 8.39e-05, 3.88e-05, 1.98e-02, 2.57e-03]    [1.03e-02, 1.17e-02, 2.13e-04, 7.70e-05, 8.39e-05, 3.88e-05, 1.98e-02, 2.57e-03]    []  \n",
      "2000      [6.80e-03, 9.67e-03, 7.74e-06, 6.78e-05, 3.37e-06, 1.14e-05, 1.46e-02, 1.43e-03]    [6.80e-03, 9.67e-03, 7.74e-06, 6.78e-05, 3.37e-06, 1.14e-05, 1.46e-02, 1.43e-03]    []  \n",
      "3000      [7.39e-03, 4.80e-03, 4.45e-06, 3.72e-05, 7.11e-06, 2.83e-05, 1.13e-02, 5.53e-04]    [7.39e-03, 4.80e-03, 4.45e-06, 3.72e-05, 7.11e-06, 2.83e-05, 1.13e-02, 5.53e-04]    []  \n",
      "4000      [3.80e-03, 4.78e-03, 1.22e-05, 1.57e-05, 4.32e-05, 2.25e-05, 7.95e-03, 2.54e-04]    [3.80e-03, 4.78e-03, 1.22e-05, 1.57e-05, 4.32e-05, 2.25e-05, 7.95e-03, 2.54e-04]    []  \n",
      "5000      [2.88e-03, 3.10e-03, 6.99e-07, 1.19e-05, 1.66e-06, 1.48e-05, 4.60e-03, 6.89e-05]    [2.88e-03, 3.10e-03, 6.99e-07, 1.19e-05, 1.66e-06, 1.48e-05, 4.60e-03, 6.89e-05]    []  \n",
      "6000      [2.19e-03, 1.86e-03, 9.54e-07, 1.34e-05, 4.43e-06, 7.81e-06, 3.03e-03, 3.09e-05]    [2.19e-03, 1.86e-03, 9.54e-07, 1.34e-05, 4.43e-06, 7.81e-06, 3.03e-03, 3.09e-05]    []  \n",
      "7000      [1.39e-03, 1.31e-03, 1.55e-06, 1.18e-05, 3.49e-06, 6.28e-06, 1.90e-03, 2.07e-05]    [1.39e-03, 1.31e-03, 1.55e-06, 1.18e-05, 3.49e-06, 6.28e-06, 1.90e-03, 2.07e-05]    []  \n",
      "8000      [1.18e-03, 1.29e-03, 1.58e-05, 8.31e-06, 2.57e-05, 5.92e-06, 1.18e-03, 2.10e-05]    [1.18e-03, 1.29e-03, 1.58e-05, 8.31e-06, 2.57e-05, 5.92e-06, 1.18e-03, 2.10e-05]    []  \n",
      "9000      [7.18e-04, 7.54e-04, 2.29e-06, 5.51e-06, 1.41e-06, 3.36e-06, 8.00e-04, 6.35e-06]    [7.18e-04, 7.54e-04, 2.29e-06, 5.51e-06, 1.41e-06, 3.36e-06, 8.00e-04, 6.35e-06]    []  \n",
      "10000     [5.98e-04, 5.54e-04, 2.52e-06, 3.60e-06, 8.98e-07, 1.62e-06, 5.89e-04, 6.80e-06]    [5.98e-04, 5.54e-04, 2.52e-06, 3.60e-06, 8.98e-07, 1.62e-06, 5.89e-04, 6.80e-06]    []  \n",
      "\n",
      "Best model at step 10000:\n",
      "  train loss: 1.76e-03\n",
      "  test loss: 1.76e-03\n",
      "  test metric: []\n",
      "\n",
      "'train' took 434.756962 s\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<deepxde.model.LossHistory at 0x7f5cf8dd25d0>,\n",
       " <deepxde.model.TrainState at 0x7f5cfbcae750>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To employ a GPU accelerated system is highly encouraged.\n",
    "\n",
    "model.compile(\"adam\", lr=1e-3, loss=\"MSE\")\n",
    "model.train(iterations=10000, display_every=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xbAdby9zDosK"
   },
   "source": [
    "L-BFGS optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 438788,
     "status": "ok",
     "timestamp": 1640284883243,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "SVSKBTp6qRX8",
    "outputId": "f7ec1668-f24a-414e-c0bd-c207bb208505"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model...\n",
      "'compile' took 0.795132 s\n",
      "\n",
      "Training model...\n",
      "\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\n",
      "10000     [5.98e-04, 5.54e-04, 2.52e-06, 3.60e-06, 8.98e-07, 1.62e-06, 5.89e-04, 6.80e-06]    [5.98e-04, 5.54e-04, 2.52e-06, 3.60e-06, 8.98e-07, 1.62e-06, 5.89e-04, 6.80e-06]    []  \n",
      "11000     [3.11e-05, 3.23e-05, 1.71e-07, 4.36e-07, 1.81e-07, 2.50e-07, 7.96e-06, 5.44e-07]                                                                                            \n",
      "12000     [7.32e-06, 1.00e-05, 9.68e-08, 2.63e-07, 1.82e-07, 1.37e-07, 4.86e-06, 3.68e-07]                                                                                            \n",
      "13000     [3.49e-06, 4.89e-06, 5.19e-08, 3.75e-07, 1.72e-07, 8.74e-08, 3.18e-06, 1.80e-07]                                                                                            \n",
      "14000     [2.45e-06, 3.06e-06, 1.79e-08, 2.69e-07, 1.80e-07, 3.77e-08, 2.35e-06, 1.47e-07]                                                                                            \n",
      "15000     [1.61e-06, 2.15e-06, 1.59e-08, 1.69e-07, 1.59e-07, 1.20e-08, 2.04e-06, 1.13e-07]                                                                                            \n",
      "16000     [1.34e-06, 1.59e-06, 8.94e-09, 1.29e-07, 1.42e-07, 6.53e-09, 1.82e-06, 8.62e-08]                                                                                            \n",
      "INFO:tensorflow:Optimization terminated with:\n",
      "  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "  Objective function value: 0.000005\n",
      "  Number of iterations: 6179\n",
      "  Number of functions evaluations: 6589\n",
      "16589     [1.18e-06, 1.40e-06, 7.96e-09, 1.25e-07, 1.25e-07, 7.69e-09, 1.74e-06, 9.01e-08]    [1.18e-06, 1.40e-06, 7.96e-09, 1.25e-07, 1.25e-07, 7.69e-09, 1.74e-06, 9.01e-08]    []  \n",
      "\n",
      "Best model at step 16589:\n",
      "  train loss: 4.67e-06\n",
      "  test loss: 4.67e-06\n",
      "  test metric: []\n",
      "\n",
      "'train' took 437.862604 s\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<deepxde.model.LossHistory at 0x7f5cf8dd25d0>,\n",
       " <deepxde.model.TrainState at 0x7f5cfbcae750>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dde.optimizers.config.set_LBFGS_options(\n",
    "    maxcor=50,\n",
    "    ftol=1.0 * np.finfo(float).eps,\n",
    "    gtol=1e-08,\n",
    "    maxiter=10000,\n",
    "    maxfun=10000,\n",
    "    maxls=50,\n",
    ")\n",
    "model.compile(\"L-BFGS\")\n",
    "model.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u9EXdIk0FjYj"
   },
   "source": [
    "Final results.  \n",
    "The reference solution and further information can be found in [this paper](https://arxiv.org/abs/1711.10561) from Raissi, Karniadakis, Perdikaris.  \n",
    "The test data can be got [here](https://github.com/maziarraissi/PINNs/blob/master/main/Data/NLS.mat)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "executionInfo": {
     "elapsed": 6595,
     "status": "ok",
     "timestamp": 1640284904382,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "aYl4TD96FrLV",
    "outputId": "30c43f0d-1147-40a7-f0c9-49d6e6466f65"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29ebQ9SVXn+9mZ5/5+P4qpGAWrmBeoyCSWCDxaKBWZlGoFFAVfozaF2KCttO20WmhwrXYARWxaLAZBWh4oCCqUIghCt2VhlYAlhfBeWUBTJVoiU0FRv989mfv9kRGRkZmRmXHGPOfe+K5114m7Y9oZGbH3jh1DiqqSkJCQkHD8kE3NQEJCQkLCNEgKICEhIeGYIimAhISEhGOKpAASEhISjimSAkhISEg4pkgKICEhIeGYIimAhIQtQET+QkT+/dR8JCT4SAog4VhCRD4uIl8WkS+KyD+JyKtF5GZbqvtpIvK/t1FXQsIQkgJIOM74DlW9GfAA4OuAn5mYn4SErSIpgIRjD1X9J+DtVIoAEXmwiFwiIp8Tkb8VkUfYtMZ6v1pErheRj4nIUwz9eSLyP710dxURFZGZX5eIfA3wMuAhZvbxOUN/rIh82JR7rYj8p40/eMKxR1IACcceInIu8BjgKhE5B3gb8AvArYH/BLxJRG4nIjcFXgI8RlVvDjwU+OAidanq3wM/DPyVqt5MVc82Ua8EnmHKvQ/wrjU8WkLCIJICSDjOeIuIXA98ErgOeC7wVOBiVb1YVUtVfQdwOfBYk6cE7iMiN1HVT6nqlWvi5RC4t4jcQlU/q6rvX1O5CQm9SAog4Tjj3xqL+xHAVwO3Be4CPMm4fz5nXDQPA+6oql8CvofKgv+UiLxNRL56Tbw8gUrJfEJE3iMiD1lTuQkJvUgKIOHYQ1XfA7waeCHVbOC1qnq293dTVf1Fk/btqvpI4I7AR4CXm2K+BJzlFXuHoSoDPFymqhcAtwfeAvzeio+VkDCKpAASEiq8GHgkcAnwHSLyKBHJReSUiDxCRM4Vka8QkQvMWsBp4ItULiGo1gK+SUTuLCK3ZHhH0T8D54rICQAROSEiTxGRW6rqIfAFr9yEhI0hKYCEBEBV/wX4HeBHgQuAnwX+hWpG8JNUYyUDfgL4R+AzwMOBZ5r87wDeAFwB/A3w1oHq3gVcCfyTiHza0L4f+LiIfIHKxfSUNT5eQkIQkj4Ik5CQkHA8kWYACQkJCccUs/Ekm4OIfBy4HiiAuaqeNyU/CQkJCccJkyoAg/NV9dPjyRISEhIS1onkAkpISEg4pph0EVhEPgZ8lmpf9G+p6kWBNBcCFwLIqYOvP3XubbfLZEJCQsKe48tXferTqnq7Nn1qBXCOql4rIrcH3gE8W1Xf25f+rHt+pd7zV39oewwmJCQkHAFc8fhf+JvQGuukLiBVvdb8Xge8GXjQlPwkJCQkHCdMpgBE5KYicnMbBr4N+NBU/CQkJCQcN0y5C+grgDeLiOXjdar6p8NZFJF0cC1hPVCVqVlISJgUkykAVb0auP9U9SckJGNit5AU8vaxC+cAEhISEpJC9rAtZbh3CiB1koTjhmQZHz9sS87tlQIQIDviY6FM+i2hhWT0bBfHSeHulQI4DjjqCi5hN5EMjxrHSeHunQI4Ti8nYT+xjxZkMjw2i11VsPulAASypAASNoRyTYI7GSn7hW0o7F1VsKMKQEReq6rfP0bbFtLgStgUcq9v7aMVnxDGmGI/zjIlZgbwtf4/IpIDX78ZdoYhaJoBJGwHG+pn65plJMQjH3iXx13R9yoAEfkZqu+i3sR8pxSqjThngM6tndvC0MtMmAZJqMVjl/rv7nAyIVZ4H0eh3/cqAFX9byLyS8ArVPUHt8hTL0Qgz8qp2dg6dt1KOW6zsqMw8PcNu9jmm1Dm2x5Jgy4gVS1F5Bu2xcwYBGV2xBVAsKMfMwG7CrYhKHbJit81bMpY2RcjYxcV1RBi1gDeLyLfoKqXbZybCORytBVAvmL/2fXZwjpR0n3WfREURwHH3VgJPf9UxsGy4z5GAXwj8BQR+QTwJap1AFXV+y1V4wrIRDmRF9uudi3YN8sAdp/nvIe+63wfFaxqrIxh542ZEfZCBsomsUy/j1EAj1qclc1AgBPZfGo2lsKUQqnUJT/7sCLL2x4AiyIpioRYLNNX+gyUjdW3RHceVQCq+gkA89nGU4tXsT5kopzKN6sAphIKmxKW1fOsz222ifaZTBBHVLu08pwA2RF3j06NTfeFVWXARmYAIvJ44EXAVwLXAXcB/p7W+YBtQFBukh9upa6NCLoVXnCxJD+rdNqwj3ckz5oU2WSKuFVvJpt3Oe7qusWu8jUdlusLq/TlZfIukifGBfQC4MHAO1X160TkfOCpC3O1BuRSctPZ6ZXKWJdgWUawxgrHWGEfx0MR/cwx/C2iiBZtI5/PZYXPqgpoV4RetuYNgas+1zpnF8dtF9Uyxlvs2Bkb22PjIUYBHKrqv4pIJiKZqr5bRF4cxd2akYlys3w1BeBjGWVQRHxGefSljMQP1TGUt6/TDHWCUOcc6nzZSIfy68pb1vPYQLD3pSyiODoW+x4eb1pGOMcK5FhhG9tuY7wuqihWUQa79K6XMTzWYegNyoMInmIUwOdE5GbA/wJ+V0Suo9oNtHXkUnLL2Q0L5VnWBVIsJUwHXlQgXyh96IWGFEIonU+LzZN56SyPvuBud9JMwm1qy/YHZfuZ/cG+qOLp4z+YbuLF53UJphhhGhKgffWHhPdQHYuU3VdWn8LII9alYhXjrsza+rCKQdhnDIbHf2BcrmEGcAFwI/AfgacAtwSeH5Fv7cgpuWX+5aXyrmZ1d+OKQMOOpQvGt2j+C/Pj2i+yEWefTernyFQ78XaXQPNZvYEY4MUKgSFrxW/bJv/SShepyIJtG0cb429dWEU4tWkN11dDwHbba1EL26+rndcX8H69MYLcF+DNOrrP2xb2wTShOgPKZuwc0KLKd5lzRUOGXghjQjhoCI4YjDGGYky/j9kF9CURuQPwIOAzwNtV9V9HS94AZlJym/yLwbiQQA6hz8oMKYCgldpK578Uv4x24zfjMkPrKocgLaAIGpa2SJ1Oqw4dEvZOEVAGOksWUDLieAgJeb99htI1aBqi9SvBMcE/1Mk3tY+87/bIGGGfRQrdkFDNRJ2StnnLkID3H1v7efPL8eu3aUNCfijOF6a23BCtTW/Xnw+0Sx4Q8CGFNTbDWNesIUbIhmfjSxiLbjx5MmfESFx5BiAi/x74eeBdVF3rN0Tk+ar6qrG860ZOwdl5vPep340zbtE3Gzbc4O18tRD304srw+X1hLPtQIVkrgz3or3B6V66K6MctEQsD6VKXUeDZsKeMG8L8RJx/M/L3NHaQjxUnl+v/VWVTof002lQOZi4Vh4/fTuPxZgCCH2kI+bedhmxeKUlLKGWyZl0b7QVj+Z+0U45fl5fcHfyijIzQtEKxxwls246V3/pGLM9KRN1wtMvr03LpewI+Qytw369XtiPq2h1uW2hXdXRP2tpKqGQEugX8uu4UaBv/MV7B9rGpG88dY0/f/x2aVmn3hCtjRgX0E8CX2etfhG5DXAJsH0FIMrZWbwLaBk3TchS99MVTtBlDTrQiLPhQmqtbYW8fdE5We2ysRYdpfPLOwWA0Nn8qllwU1q7kxxq7sJzzV39vkAHOCz9dHUZ7XTzMqMosw5NtQ5DJbBdOquASnHp7G+hXZoqlLYcM4ZVxYXx0neGuJ/O0bzggrOCoLUvfrwNqCO7PKIu3tJEPAWRqYvrCPGsdMrIXoCYZ7X6nHk0K+xtuTOPNssKV26bdiACme1F1XueUXQOFOWUHNhyTGMeSOGE6IHYOkondA9k7uq1gtpPXwv+WjnUisRTEG3lMRDnl9fmvw9DSsLHmCAdkht+/j7D0qYJypeWd6AS7E1Z0pA5ngwY27QSowD+Fbje+/96Q1sZIvJo4Nepet8rVPUXh9LPKDk7OwPEu3w6fmjp5mta8VLTpBlfCfGmNW3pUAv7QjMOJXfhqt6uUjjUmevQuSnvDDPnxgnNVFyd1Na0LfdQcyeAD42wn2vOYWnDJq7MXdimnwdoh2XuhLilFSrMC/NspRXsnlIwtLIUL2yEeSmooWF+tRAXduNUxYXFi3NhK19LP1zHefK3kb6PNopWl1HxaD5LWePRIFM3YXPjPlOXTjPt0MiNUsgUyY1gzz2lYPLkJi7PS6cMZoZ2kBfOUj4wV6ccZIW7RsWepp9lJSe0Cp+0iiDzZyG1ieEL/qq8eR32fk8YwV/T5h2lkEvJgSnbF+Iu3psp2HBTyDdpTddSaA1iYOE6QBubG4RkT9N6b7t/pTNDLZCOS6eiZQ2aX9+hzup0LUVR+J4DE3dG80EZAnEK4CrgfSLyh1T9/ALgChH5CQBV/dWIMjowH5Z5KfBI4BrgMhH5I1X9cF+eDDhL7HGM8RFcQudyqkJb8ZaONH6r+CatpOykq5SCtVyqxs6prR77Qg/J6x035kUic9wrsL57SkpjjVlrptD6ULlTIpo1BD/A6XLmLPZD2wnKmVMAZ6wiKDPOlFW9h0WtHGz40Ar9ImdetBRAkTGfW2FvBHshqEmnc6s1pfoDxPu1N3nUtOoPIJsHaAUdmv+bmRfq05wr3P1ql6a46cWoUrBCvjbxXdgX8GpMdvu6ylxcuP7NME3vfnUGpY0/sHHqukk5MwvxM5gfqMljHmSmyIHpf4Y2mxXMTPjEbG5+awVw0tBO5nNO5eZdmhP2pUrt2jGKJctq6/zAvMBTMncC+1R2aGiHDcEPcEIKF7ZxJyhcuFYOtUvphBtP9VqW049A3jLiMiBva2QgCxh7y6L0ppVFY6ebhd0soZ2ZeUktd2z6wnOvujU8T1EUDZo1HA8NLXNj3jcwzxhaaTwNmZbB2YiPGAXwD+bP4g/N780j8g7hQcBVqno1gIi8nkq59CoAEeFAhINAXJ/WLlr+AN/H6+exu2YaWxXtgqv3ouw++Ny9KPVoVd4zQsNXD0Yr27CtWTMXLgI3h9Q++azrY0dqK99ztfiCHyor3gr+M8XMxNWzgtNzmy5zlv2ZubXwa2FfGFpZCKWhYX8LQYzwdkL8UGrhbfxX2VzIXNj8HuJoMldHswe+s0PzXg7Vhd3vmZL80FjJ5jc7LJHDwtAMA/MCOTQVFoZWlC6sZVnHtf1HpdadxvlzMsQITnLz3ma5C+vBzP3qSdNuJwztREZxYN7XSdO2JzOKE8aFdtK83xNQmLCNK05W9Cre5D0B5QkjfE6U7vfwRPVsZw6q+g8OCk4eVG3gFPws61immZTMTB+aBXzOvnXuC36oFIENO+XgKYUTWBdQyYFVKG4dATeurYD3BbsV5jlSG1ouXS3khs6p5BJwzeiwvV+6KaY3bl1e7dJEnbJoKAqp89g4u4BvlcOhPyvwDE1bxxn73KpOz9n3UXh8+OEzg08Xtwvov46lWRLnAJ/0/r+G6ubRBkTkQuBCgHPP6b9eKWN86gaVBWE1tO0ODUUQoC2K3Hjyt42+bZb+om6bZqEqblHUGcTq+eW9dLX5S/1rG8y3tFtunIZ13rDsjVCxSmEOmVUG3q8v+AHywxIx4exMlVkOi1rwH1raHM5UgknnVinMawXglEKBWgUwsEIseY4a4SNG2DObIbOZoR24OC1MeG6UfzGrTX8PKlaoGYKIcxH5MwtnV+R++5k81q3mudrUraWUnbWWEGLPzYztsgnv1um6bFxcIJwjDcFfxWVRgj8k7IN8jiiFzAlk9dbm7JqGuBfiKwI38zCPWUR4K8C0mfV4jqxV2fYfkzNjaxwxM4BJoaoXYT5Bef/7HehhZ4Wvi9DiaMPd0zbyGnmtXz3sCmqsEThac5HmUGfBRZ+hHUQhZG4qXM88Mm8boOsEpsPNssL53e2CX5F1WyMjsOgoSp4ZIWQEzzzLnM+5MC6eIs8oZrXrB6CcZ+iBES7OPYSbFfizg7YbJ5vXs4eaVruInMto7rl7rKIoagXhu4ysQrHKJpsrYgW65wKqaS1XUA9cl8ikFsrOBSS1yya3Aluca6ecmbaa4bl2zO+BUB60aDPPHWTcPuWstvbVuYIUjAvIuoJms5LZQWHCxvqezTllZgDOFZTPOWWmWidMPzmZzd0sOCTk/T5uZ81uYbbESW+b7lTAr30ghZst25lAplq7frzZuF0msaZfLuLRmrMDn+bLvWwJY6xt2VfP1OwgIbdQU85Y69xz/XjdsC1rysa6QC1f6jY3M3OyQZnj1gq0dhX1YUoFcC1wJ+//cw2tFyVwY88g7XUBjSzYVHnr/9svwKeFhH3oZVT7b/tpzlfnreaHNPmQNZZLvTMjNIXxtwnOW7s/5mXOCesqMm6heZ41Fn/BuIDai8Bl5tYFbDsWRebWA4qitjxLf6EXmJf1ugD+YnB7wbeQjs+eUpr+e/PcbnbhLSDXYS+9DtDG0F4DgHoh13M918pAXRptp8tAc2vZm7i8dBLOxpFpvSDsfksXtv752ax0C8L+YvDMCXnjk88L11+ssD+RzxuCvyqjYOZ29XgCzm0gmJn61b03a4AcSuH69glTxo1yorH4C3BAdwdRjr+DqF4fiFnwbWwHBS+eUSy1CNyzjghNeVO03LZ+fLVDRzp5Y+RLyMD0dwH5imLsuokpFcBlwD1F5G5Ugv/JwPcNZSgQrtcwy+MnfbsvoZm/a52HDkuFzgsMCXt/C1f7nECh0kl3qPng1i23f1mzemeGpwjaB3VmZLUCsJ0lK9yOH2sB+ruA3HbQwM4gVXE0f5unUwY+rbULKLjn39sa6hSGeq4L++Cei8oStZSWVDYu/LaLqh32yuzFyPbPes3RCGfx8pg4EfXCNc2+I/G3gZpw/VuLjXobqDa2hEIlBJ0ycDRvZ5Cn/Nu0g6xohKu4Ini2wV97ArixPKDImpbpgRRkcuDCYLaQOgVQ98322YBc6nQWOVrvkhvY8ukrgLEtoaF0sRiasYd27dT/Z4PpQltDx2RJuzxf2Ne7hbprPG30KgAR+Q0GbCRV/dHBkkegqnMReRbwdir751WqeuVQnoKMz5XhTxKM+S6HXl7voa+AsG+nKwIvtPmS+4V96JRuSOP3PZsdFKXaXUPqdhO5mYB6AgIrzHMO2tPPrL41NHQ2YOiA19w7a+D/6kBe/9BX+7CXn5dBWtefHfJvr+sqiKFDX37Y32jQPszlKwCXxov301mBbRV98LyAaMdi9w+M+QfC8gFaPWPsHsgqyGpT15thWoVt8x5KHjz01T4AVimArvXePgcQOvTVd2I4fJXEgALY4EGw4EHTiAOkFb2VzpMbflnt8z59Nw2ssgvo8sGca4CqXgxcHJu+0IzPFDdbqI6YezvGZgRj10QMndqry5NOXMga6D3RF5pODrmInF/FS+MpB8vPzHpitD42Xri4epCETgfX/B8GFUUnr3Y7fejah7GrHhongIPvbj0CfwjBax8C9lL7IFnvFQ+BMtwaUCidJ7i7102EBXu73FyU9lUQpUrdZ/wdMvYQo5mZNYS4m5XW9fmngy3ap4Srcvot9r640NpEzB1EIayiCMZkS+ytvmNXRcSc8A1dD7HSbaCq+prR3FvGXDM+V5zVG7/oJU2LXru86p0eUa4nDd3J4888+hVA6ORf6DlySjeg3S2eUjqntB1MpagT9v5VzbXS8BREh+dwG7T5WeXit/FbFhdf/FsUQ26G3vt3Yu4MGthF0663ffipqUSGXSKu/RttZXeY5K68vLVdshRxU7YsYKn38V3x6/vsA4J7XRfDTfyFtKH+F3vb8NhFcUMXxMVcXR9zF9DtgJ8C7o33SUhV/ebR0teMgpzPzBebAQyXt7oAWeYmv4XTjVoLwxZzjFIoe5RM+zZQ/16ixsBujc9M1JUzdL10u26g4Wv3bzoNYdjKifuC05AiWflDKpG+5qF6gtcbRPJVK3gN9Gd/77h9V+pciq5+VeatcyqxVz/3WeHDzzv8bDGW/SLlrQOLzjoXueZ5KG5cHizvArL4XeANwOOAHwb+HfAvEfnWjkKFzxc32XAdcVbjKtdLj5UR+6JjeFnkjvBF7+gfrDdi+mkRtExDSiaENYztqe6Tj623bgtPYA8uYIfKqInu3qnAFsF5JF9j1vW6P0QTW+8uIvrrXiNjZvhK9sUVCsQpgNuo6itF5MdU9T3Ae0Tksoh8a0ehGV+cn1x/uSsc2lpE86+rI8CIAFigzlE3SqQgX+9nLP306/HnT/2RGIvGNd4LPltIYENXYIfOsPXyExKoEfljD1r11rsOzb1HWLb/jY2rVcd31Cchze+nRORxwD8Ct47It3aUCF8s1q8Aeutbow95FUG2yc/NuTqWeNZtP9O6eZgCG1FEAV98fN7l+/hxE+KLYh3vejEDc/H6YhTAL4jILYHnAL8B3AL48YVrWgNKFb5chG4CWmcdm184rOtarzDYVaG6DSG9b4pgX7Hrn19cB7bVl9Yha1Yd8zF3Ab3VBD8PnL9SbSuiVOHGDSuAZTC18Nlm/dtUkI16d8SFc+yx5/J/6rFqsQk+NjIDEJF7Ab8JfIWq3kdE7gc8XlV/YXEWV4Mi7uqCKbCvQmhXOn0Iu8zbpnGcn33fsKnPi45h0zInRpq+nOqrYL8FoKpXiMjrgO0rABV3tfE+4ygP/KkGyi5hXw2Fo4Z9GWdT8hmjAM5S1b+W5scV5n2JNwmlvqgsYXrsywDbNaR22xz23QDZdt+IUQCfFpF7YLx/IvJE4FMb5aoHirhLyRL2B/s+KPcRScksjn1vs2WWZ2IUwH+guo//q0XkWuBjwFOWqGtlqNZ3kSQkrIJ9H+z7jj1fS14rdtoFZD7Z+K0iclOquwBvoLq6+RMb5i2INHCnRRq4RwtpPK0X+zbbHboO+hZU1v85VN8Bfqf5/znAFVRXRGwdQwIodeaEo459EzBHGUdB3gzNAF4LfBb4K+DpwM9R3TLynar6wS3w1oEiyQWUcGxwFATMccQ+KekhBXB3Vb0vgIi8gmrh986qeuNWOAtB06BISIjBPgmhhMWwyF1PYxhSAPYOIFS1EJFrJhX+jpfUsRMS9hnrFGAJq2FIAdxfRL5gwgLcxPwvgKrqLTbOXcJOIw3khIQ47KrhOvRFsJ07caUkoZOQsAnsqoBK2Cymu1hnA0idOCEhIaHGmEzcOwWQhHxCQkIMkqwYx94pgISEXUISMgn7jD1TAJIGXEJCQsKaMMmpKhF5nohcKyIfNH+PnYKPhISEhOOMKWcAv6aqL5yw/oSEhIRjjXSvQkJCQsIxhahuf2O9iDwPeBrwBeBy4Dmq+tmetBcCF5p/7wN8aAssrgO3BT49NRMR2Bc+YX943Rc+IfG6Cewin3dR1du1iRtTACLyTuAOgaifAy6laiAFXgDcUVV/MKLMy1X1vLUyuiHsC6/7wifsD6/7wickXjeBfeETNrgGoKrfGpNORF4OvHVTfCQkJCQkhDHVLqA7ev9+J/vj1klISEg4MphqF9Avi8gDqFxAHweeEZnvoo1xtH7sC6/7wifsD6/7wickXjeBfeFzmkXghISEhITpkbaBJiQkJBxTJAWQkJCQcEyxkwpARB4tIh8VkatE5KcD8SdF5A0m/n0ictftcxnF50+IyIdF5AoR+XMRucsUfBpeBnn10j1BRFREJtvGFsOriHy3adsrReR12+bR8DD2/u8sIu8WkQ+YPjDJlSci8ioRuU5EgpstpMJLzHNcISIP3DaPHi9jvD7F8Ph3InKJiNx/2zx6vAzy6qX7BhGZi8gTt8VbNFR1p/6AHPgH4O7ACeBvgXu30vwI8DITfjLwhh3l83zgLBN+5hR8xvJq0t0ceC/VOY3zdpVX4J7AB4Bbmf9vv6N8XgQ804TvDXx8ojb9JuCBwId64h8L/AnV1/4eDLxvCj4jeX2o994fs8u8ev3kXcDFwBOn4rXvb3QGICK/FENbBiLycaPJPygilxvyg4CrVPVqVT0DvB64oJX1AuA1JvxG4FtEZNvXhI7yqarvVtUbzL+XAudumUeLmDaF6lDeLwFTfvs5htenAy9Vc3pcVa/bMo8Qx6cC9tOptwT+cYv81Uyovhf4zECSC4Df0QqXAme3tmpvDWO8quolWt8aMOWYimlXgGcDbwKm6KOjiHEBPTJAe8waeThfVR+g9cm5c4BPevHXGJoPl0ZV58DngduskacYxPDp44eorKwpMMqrmfbfSVXftk3GAohp13sB9xKRvxSRS0Xk0VvjrkYMn88Dnioi11BZgM/eDmsLY9G+vCuYckyNQkTOoTrn9JtT89KH3nMAIvJMKlfLPUTkCi/q5sBfbpqxowQReSpwHvDwqXkJQUQy4Fep7mfaB8yo3ECPoLIA3ysi91XVz03KVRffC7xaVV8kIg8BXisi91HVcmrG9h0icj6VAnjY1LwM4MXAT6lquX0HRRx6zwGIyC2BWwH/DfAXuK5X1bFpT1zlIh8DPks1Vf4tVb3IDJTnqeqjTJo/oJpu/5OcPPH1B3fo3GeUkJCQsH/Yok448/FrP62By+B6ZwCq+nkR+SLwdar6iQ3x9TBVvVZEbg+8Q0Q+AlwC3FNE7gZcC9wDeJSqXnnyrufqHX5+V2fRCaPYTSMoIWH72PJY+D8/8NNBGT64BqCqBfBREbnzJphS1WvN73XAm4EHGZ/+s4C3A38P/J6qXikiz98EDwlbgJCEf8LxgUT87Qhi7gK6FXCliPw18CVLVNXHr1KxiNwUyFT1ehP+NuD5puyLqRbNHFT150/e9dz/skqdCVvGDnX0hISFcEz6bowC2JTQ/QrgzWZxZAa8TlX/dEN1JWwLx2TgJOwoUv9bCKMKQFXfs4mKVfVqYLJTfAlrQBpsCasi9aFJEXMQ7MEicpmIfFFEzohIISJf2AZzCTuAPfBjJmwBMX7tZf4SJkWMC+i/U1238PtUe9n/b6qDOAlHAWkQHj2kd5oAIONX/Ud9EEZVrxKR3OwK+m0R+QDwMyuyl7ANJGGwu0jvZj8RIVj3BTEK4AYROQF8UER+GfgUO3qL6LFDEiDbQWrn6XCEhO0uIkYBfD+VwH8W8OPAnYAnbJKpY48kcFZDar/1IAnftWPXboSI2QX0CR9NvzAAACAASURBVDMDuCvwB8BHze2HCctixzrB5EjtMYwkiEexa4J1XzCqAETkccDLqO4+F+BuIvIMVd3ZW/gmx3HpjMflOdtIAvl4C9wj9P5jXEAvorqy+SoAEbkH8DZ2+BrWrWDfB8C+89+HIzQ4QzjSgvcov7sdfW8xCuB6K/wNrgau3xA/u4UdfWkd7AufcGQG+ZERxPv8Pvb4HciOtHuMArhcRC4Gfo/q2uYnAZeJyHcBqOofbJC/7WFXOtOu8GGxIx11DHstkPekjYHd658t7Ipg7WBH2y1GAZwC/pn6Yyb/AtwE+A4qhbC/CmCbL2WqDrADA2LnhfMOtFEUdrQdd0LoTtw2O9HHl3gPMbuAfmApZnYVm3hRm375Wxpgk3fiXRAkQ5i6fXqwEwIYtt4+W+mvW2zbKbpX1EngvcY6WnVdb2bNnemoDYBw/dNWH4udEcIWE7TbRvvjhtp33SyvvR9suF8dTQWwzFtdpSes8JLWOmg20Vkmn1rvmGBdBDukvPbRWFiV5bX0nYnH9qb7f8w5AHsH0O4jtsEXfTELvoSFXvwqL3iNg3rrgnaHhOO2MbmrzceOCW1Ysi8u+RyLvItlx8h26lguX8wM4P8TkTcBv62qH16qlm1grJFH48cbcPRFxr6EyA6xklBek5DZKWGVMIwtKPGVJsqL8rdA+ph+Glt/bJ8fK29d9cWXs4FFYKqPtjwZeIWIZMCrgNer6m58E2Co8UJxPY00+BKGGrYnX9TLGH3x40U0M6wwXV06535h5x1KGit9dv5JHKIF0xqMsKG6+vIO5xkT8svy0iOHBurKVuCzDzG7gK4HXg68XEQeDrwO+DUReSPwgtYhse0hKNxDtGbDBF9YqPEC6YKNHEwXqKOvnv5ixuteoPzRrBNqAF2nLBsRoKFYjRW6fp6Fc/iZ19DY6yhjEZi+FfvcPneh9g3253a6hWYA3bSLCufYMqLzeuGQ8F5UQWSRvAwpijai1gCAxwE/QHUh3IuA3wX+DdWH23fj4zB+43gN0HmBfuN4cZ2GbMSF6tNQ0oHylu9cvTxE5l023apYVLDGph9SGH1xQ3mGal1GOYwJ52iFt1NGfls4L5a8GadhpdAmqHT6qio9iqGZW0Q77Rw7AxgblzIw9rO+dK06xoR5Xzmh+BAfoTQhRK0BAO8GfkVVL/HobxSRb4rIv1nYJw8J/YCwlyDNL6/7cpt5+juLTReiNfP0K6hlrAsfi2j/sTp8LCMIy4E8Q7WF6grTQvFh4dIWBn4a9bnRkCBp0hpFjdUVMnRD7RJqkKH2W0U5hMpdot8MGSWutOCsOTwjqPP49MCswPAvbbpff8CAU/X7+aDajzLOQhZ+SOhnAXnQJ+BdOYGyQ/+HymuER6yIQQVgrP9Xq+rzQ/Gq+qODpW8KQpzgF69xfWEfaOQ6XS3EhwT6WFyMZdB88QFeBvL2xbexLUdBiIOQAmgLPz+NH1cGBHYob7veUDrVumz727QQxaVzfcdL33k2lTpvSAlpM60ru0VrpuvW0RvXS9uuW2hYyAcy9Bha0GJ9KaNOWrQ67fg4b7EpEhij4sVXNKXnMVuCPyjgRTvCPiTE/XQhwR4qO5S+D4MKQFULEfl2IKgAdgmSadfKl27HyUQbL9+ml1beEC3zX1rjBbXqpWsRZJ5SCL2gUB2uLMLp2nU18oRoa/IplIFuHxL2bVoZEuKISxdSBn6ceuXYNGVLsBeldPIqXQVQlr4VUQ9yLZvpFGohXnq8OQVg43yaFxei0aaBtJWC0hXyKl1Z26NEogz6sTRDFv6ogA+kC7lfvV8NuWHbulK0k1cDNLzxO2gQZvW49JWDnXn4hl6WlSbsxZl429can0pcQh60x3cm2hH2voxwvyiZlN28Ix0hxgX0lyLy34E3AF+yRFV9f0TezaFtgfsWe1anybyXCpBldceoX2hXsGaiLt4X8HnWbGQBZq6cbsPn3ktZ9KXZF1/RSlNe1wqo85YuLpTOh592WZRad/e2Uiga1nxmfsWlK5yQzmohTy2452XWzItPs2VkjmYF9rzMXNmFiStKobRhx1dG6Qt0wjMFLamVgvmllOqvYswUIi4sfpydcbg4XF73ikpPdtp0vgLwu7C28mornmaeBmL1/5hwpy3Y69+OYhCtac6i8tJlHmsdwQ5k2qxPQL3x7dLbciwt8x43Uy9vVwE4WeLF2bCrNispy7wqzsaJNsKuOvNbmHeZZ3U3ySU0FusyrCzxhb6N9+PacmOWlUF5EKrPR4wCeID59WcBCnxzRN5BiMijgV8HcuAVqvqLcRn9sGexZ7hwTasFf/Vbdl5a7gn7PKsb0Qp724h5VtYNbgVygJZJWb8sL86+tAPv5R1k1Rk7+/JmWUFOVykcSGF4qTuBH7ZxTlFQdwKLMK2rCBZdR/At98KzgazwLpxgzxphgEPNKbywLc+G7e+8zDk06eZmIJ4uZ8wN7UxR0c6UM6cUDk26wyLn0NDmJt28UA6pwo31AKuMCiP0C4G5eSZDk0LAHI0URwOZt4R9AZkXXxXcUgY2b5vmpauVgiLWN9ZQCnRpNGmj8CcmbWvaCzcEsaU5C8lL5wSyeMLZiwvRWnkrgS3N8jKt89o4XwFYhZF1aWTNMJjmcXnqOKcMrKzIxNG0rBVGWTbli2aKWrlhHqMoKyUAXl8bmbX7xmRbKcyyIiBzurLEp/UhZhvo+WNploFZX3gp8EjgGqorpv8o9rCZtDup0FAGUL0gX/AD5HmtFX2h3xb2s6x0NF/YzxytcOls/Mw09oHUSsG+gJkUTojbMg6kcELcxh146TIvzikALG3ueD2QuYlTr7yalmGfrVYOtaLo0ixyT3rkkZKk8KRGW9gXiFMQPu1Qq254xgn9mRP8N+qBiZtxujwwtCr96fKAG4oTANxQVr9nyhlfLqp0X5pXtBuLA748r2in51XeG2VWzzychZ/Xgv/QjNjDDDk00/szRpgfClnVvI6WHVLTDvFo2qTNPZpJL4V6NNPOJWSFungAKT0F4GYe6s0GvHfUfl0NYe4pvJZ1rrl4grgWsJpLg1YJU5unTq8t4VzmnvDOa2Fe03C0Mq/Ldnm9eFtGXQc1L510dHjRDDRvKYBc3RRfPaFfP6eh5eopF/M+1Hex+GPHFO4pglJtHnG/Y9Y5NF1AVr5koh2Zk3s039Acm+lH3QVkPgv5tVRXQwPQtzC8AB4EXKWqV5s6Xg9cACx12ljE6+O+f95qT29K5yz/rLbsfcFvaU54e4pg1qKdyObOovcVQB32BLuhnTTC+SCbN+Lt7wkb72hzJ4CtYD8hhQv7ysGGT3iziAOayu0ApTXWOCHiwrlpyUzEhS0yz8LPIw8PFOYmkZKSwjyHHRAFyqEJGxnJocKNThlU9d2oM6cMvqSVYP9SeZIbypMAXF/cxPye4otFRfvizMTNT3F9Zmhy0tQvzAszU7AzAYzFD2AUgJzOyG+saPlp89w3Cpn5KvbsRkM7rXX8adOXTpdkJpyfNu/gdEF2pmoPOTOvfw9NeG4stsM5OjcaojC0onRhtcK+HB7gZOY5RGpfZm6kpYgLi0nHbFabq4amsxxmpo1cugyd2Xjzm2fozAp+o+gPpBb85reciRPYpaPh0TDp6aTzaf5v6SmhbjqtaUVTaVHWYfGVUmtWpV7YR2mHQ9kdC3atoKQW/Pa1qYpzV5WeUoiZffetE3Y9DLU86EPMOYCXAWcB5wOvAJ4I/PUol+M4B/ik9/81wDcG6r8QuBAgv83ZAf76G2yVk399L6Ltnw+nKYNunNhF2CzglnHWu3P7lE4p1LOCwnPzWGFfNgR/VQYcuHKNsKcW/AdmJISEvU/LIvcY1W2ZUxpeD/Gnq80pczXSmh23oHAzigMzYk9IwaGnJMEoVTNDyErf/db2j2q4H7QWY6X0XM1lHWdn1s46L+hY8dmhkp8x78sI/exM0RT8gJw+dAqAw0oN6nwORgHoYa0IasFvJUnpXBIW0thSZoRzJogT/KbePIe84k+NwBbVWpoaoS/UCkdmtbKsDS47U9DaNWbbeQ6l8/3bmUydx7ltSy/elpfhubXs8wZ2YCnNNRH764StR2t1Wd+Fpo4mza3BfmQ77GjSiRxy96wLMbOIIcTMAB6qqvcTkStU9b+KyIvY4veAVfUi4CKAk3c9t/O09dY7rRd2qEna1ryibkGwnrZlbrqWlbWAI5CMomqyMrNWrVAaoVx6Lo+5nSkY63cuGYdZNXicZa950AV0Wg5MuBqomWhnVuBb+77bpz1T8F07J6zQFV951AIxpzmKMubBDjbkDmq4gAI7fmz8YcMFVLt+oHIFWav8RuPaOdS8ngEYq/90eeBcP74rqHYBVelumB9wo6HdcGhcRUXu1gMsBBArEK0le7L0hGiVLjuA8oRRiKfM77x2B+VnrALInOsnP6wFYmbDjqbIvHRhqBSLmFVE6/aRQp3FbxUP6m05HjscUG9rq57N8887IZ7VMwX1f/MuzVn03nSynLXS5U13UBXnzQYaNFye6lc6rqIyp+EisnnppNMOTb01gIbLyM4Q/PWB1hoAeb2e6H6lNih8WnsLqW9shLaKL7Pm1t5MgXgbJmx5mlGsYRH4y+b3BhH5SuBfgTsuxHEY1wJ38v4/19AWQ0PjW1+u3TkSmh5nbuqVOeVQdnZ/5CocGgHhu4Wsn01k5uLa/rhlfHS1318b4SpvLexDC761Uuj69v3OlXtrAe1ZRljQx+0U8hd+G7t/sDtvMpfO7eChXvj1F4Ttr81zupw5mr8gbONOl/XiL1SC3YYP3cJwXofN7+E872whFVGymZmWG2WpWVYJE2qhL4XUC77mNyvqhV4X59H8xWCXzl8EdrTaSg4tDPszEzCWcevVNRaGY2HtKAnRpBPfWLR1kqwVjxH+HRreAm392/HtS4DmC/FQeX5doXQhYe/tHKp+1Vv8rX87m0y8DSWhnUHNvfw0aCGE4koVN4PyN1tYmh0L1ZTUakZbXunWt/oQowDeKiJnA78CvJ+qa70iIt8YLgPuKSJ3oxL8Twa+Lyqn1sZOfcpPPAvIEMvM+ZrtFLYstbFiD1XDFnY137o/stK9tLm3zbMWyp7l3NrxkxGg9SzmxGz5rPx7Nr6Oa7uUMs/VYeEL8bGzA4tiaO9/czdQPTOK2QZ6WNbC2e7yqbaG5h2a3eppaUVZbw31t4MWZiBYq79U6Q4Oz1LLzMgopaytX7tYXIg7E4CrA7cjyG3zbOzk8Wi26Z0Qb4Ux7qaWsA9u+Qy8g7FdQEOvfmx//1C8ipd2lNYUzo2toYG8/tZPP4+Lawtx8YS81HnxynHp21a8f6bIXztsrycKja3kVVy9ITr35MbQ4U6L0lsDcLvRqMfP3MyKM+2O47IUJw8KZ3Rmq50ErhjRF5jgm0TkrcApVf38WL6Icuci8izg7VTu31ep6pXxBdT+WBtw+7Tdy60PD7nteqINZVDRAtO2LKtptuEDU7lMoL1u0HfKL3QYJLSmMHTYa+iYt7/iH3s4bBWMHf7y9/C345vpmpZ46HCYf+jLP+BVtM4GVOlMuV6ck7nu3EAzDDROj1pkuSdIbP/K8WabdT9049IdJqMrMUuGD325BUFtxht0fN1enk6adrp1IKAAgoe9QkojOLvw8rXTCfW+/mAZXpyM08YOizbGtgtXcdX5IW2l647VPhnRHvs+/ENkblzYehVn0Ph1WOOlNhxr488aQzHrjrG7gB4K3NWmFxFU9Xdi8g5BVS+mulBuhTKqX6mWqiqamwGo93Zt+nqm4PpH4IX3nQ7G5fHTNcuD7gsPXfsQe/lTI66HPkRrY+zenyGM3QkUEuyN/APp2yd42/W5dc+WwvDzaE+8BvK2aUDdD3wm3fvyyrJ9yPqSVWgL8b7rHzybpRnn0YJ5+xB6nesW/G2MzBTCeXxNFsjb7pcSShc44euVF7wFwM9LO11onHtZvPEbvC3ApMsa6WrB36iLpmXvrHf//VqjxVcszp6tjcWukZh3NzpEdIKYXUCvBe4BfJB6a4YCKyuAleCezQ58ra2n0BTSppZaKdQ0CQ789r0hftl+mqEL3UJKYzQdXYQEe4wgX0XYL4MhBdEXF6MoQoe1mvH99YQUQThdlzeRul+FrNv6YjL/biFtJO/UFRTYw880lK6bZjzJwoiodqFvbYSEvUvfryjGLm5slNGqI/bixj5Dr5OOYViDJUM6q2kZdGhKdyyEvAB9huOQN6EPMTOA84B7q45tMdgilLr1fUXg/GeW4lllvm8tYBk4RdLoEO3eF+5wwUYe6HxD+Ra7IjrulUz9Va/YnjM2uxhWLqunj71MrZ5rEjYOWtP4KtjsmzEMDXGz1DXVKyDakBhga9GPLgX1RkM5DI+fNs99ccumCyH0XkrPOLWzgkK7V10XdA29Xi9Bq45lvQExCuBDwB2AT0WknQ4K7WYJXfMr/uj1lYJFwFCrX0KtPPym7SgKH5Eduy6r56VtWNhva6awDsG1sCmyQJ1rEawDbSl9b3/BagN3lA62y0aMgNg+ORQ3VMYCM4pVjSkfWjsJho0Ir7zCzQTp9DdfEFurXzyj0+fNbpNex8dfYp41RgHcFviwiPw1cNoSVfXxEXk3B08m99KCroUepdAuA0YVRKfaYIN3rcFgqogO5/MUi3WN+227kiy2belu/CnH2nGF5w36xteMZbiLnz3EuDTj623vFKxoXaEb7mMhBduVF6HyxksCtNcUAJYT7MuM0RgF8LyFS90mhhRBmw7jSqGRdUDau0Q1YVBYjXTcYSUSKChyIAwpm0Wwact4UmxZyawFW27LrbRQe60llMT3AA0I4laugfIkOB6GBHvvWlEvH+E6+sqA8LpYKN1geRFpYraBvieqtqnRnVHV9DYilEKVNaKhezrQIB+DHTzSRbDA+G93mLWv5sRKh30UtD2Yel1lpbZcQnkM5RjjJMZ4aPTRyGeLGp9e/X3vbEigt8vow6AbOJg+RjD1pR3jJZ6PXgUgIv9bVR8mItfTtalVVW8RX82WEJoNDKXzEVQeMWq7++3R0SzL2FMDM4/RrFHPsRg7zQrik07lSloLGhbndGyMYVQArHk2t5amWIYnL0+sbz+IkNsmdkOAc91GJfdzLjAWurysywDpVQCq+jDze/P1VLVFxCqCUB4fUZb94j77VYRHZ6fJuqBsxa3QGVhTW9IjaFqm0/HRi6D7YoPVbaj/VXs4lmc86iPzvZljhf0ax3Lv2sOSxS1ZVMw5gFsHyNer6mGAvlsYWgtYNH8bm3R7RPpAV0Gww2zTRWOfcdeEaqsJtr0IPYTw7aXb5WGj1XX87isWF5No3TNpv+jBwtbTr+q1h+XyxywCv5/q0rbPUj3T2cA/icg/A09X1b9ZruotY1VlMFReCCu5VJbMvFBnXq6KURZ2fT1gdBfOdtgYRcTC5DYhEyhsXfeYDVayPn96p+jls0aP5VX7RIwCeAfwRlV9O4CIfBvwBOC3gf9B4A7/ncdiayrrqWMMa1FKa3qQlabi62FhCCtZhruqeNqYUhHtmPIBgucetoHB/rzJJlnHduAIxCiAB6vq0+0/qvpnIvJCVX2GiPm80r5jWf//pnkYwo52viDWvLaw7QXYtSy4Te1KWuQd7OAsaGoFFERjS+quNNpi4yNGAXxKRH4KeL35/3uAfzbf9I27MH4f0btjYKtc9GOZ/jYV75savFs7vbyVahpY+zbTXRCg+zQLisGuK6gIxCiA7wOeC7zF/P+XhpYD370hvnYXU00J14FVBtQuPts2Bt1kJ6AnqRbY4BmHXRWSy77jXVdQEYg5CPZp4Nk90Vetl509xz4rhzGso7PvYxscYSXTh10647CVA3e7qpgWwZJ9KGYb6O2A/wx8LXDK0lX1m5eq8bhik7uG9gXrEixHra22LYB2TOEMYZeUUQiTnwi3WLIPxbiAfhd4A/DtwA8D/w74l6VqS+hHUhDxSIpkNUxt8e6RAhrDriuoMcQogNuo6itF5MfMvUDvEZHLNs1YQgsxHe24CrRlMdWdSMcdUysgH0dIGS2DGAVgT/x+SkQeB/wjEDodnDA10ixiWmxKlqT3tjnskjLysSXFFKMAfkFEbgk8B/gN4BbAj2+Uq4TNICmI/URSLMcPW1JMMbuA3mqCnwfO3yw7CZPiKO9iSugiucCOPWJ2Ad2NahvoXf30k38RLGG7SMohYQybmKmkvrVRxLiA3gK8EvhjjvLJ34TlsQtXaSQcTaxTqaQ+2UGMArhRVV+yzkpF5HnA06m3k/6sql68zjoSJkZSCgm7hqN2En4NiFEAvy4izwX+jOZH4d+/Yt2/pqovXLGMhH3CNq73TUjYBBZRHnvUt2MUwH2B7we+mdoFpOb/hISEhAQfe3RmR3TkKJuIXAXcW1XPrK3SygX0NOALwOXAc1T1sz1pLwQuNP/eB/jQuvjYMG4LfHpqJiKwL3zC/vC6L3xC4nUT2EU+76Kqt2sTYxTAW4ALVfW6RWoTkXcCdwhE/RxwKVUDKfAC4I6q+oMRZV6uquctwsdU2Bde94VP2B9e94VPSLxuAvvCJ8S5gM4GPmKuf/DXAAa3garqt8YwICIvB946mjAhISEhYa2IUQDPXXelInJHVf2U+fc72R+3TkJCQsKRQcxJ4PdsoN5fFpEHULmAPg48IzLfRRvgZVPYF173hU/YH173hU9IvG4C+8Jn/xqAiFxP/25uVdVbbJKxhISEhITNYnQROCEhISHhaCKbmoGEhISEhGmwkwpARB4tIh8VkatE5KcD8SdF5A0m/n0ictftcxnF50+IyIdF5AoR+XMRucsUfBpeBnn10j1BRFREJtvGFsOriHy3adsrReR12+bR8DD2/u8sIu8WkQ+YPvDYifh8lYhcJyLBzRZS4SXmOa4QkQdum0ePlzFen2J4/DsRuURE7r9tHj1eBnn10n2DiMxF5Inb4i0aqrpTf0AO/ANwd+AE8LdUB9H8ND8CvMyEnwy8YUf5PB84y4SfOQWfsbyadDcH3kt1TuO8XeUVuCfwAeBW5v/b7yifFwHPNOF7Ax+fqE2/CXgg8KGe+McCf0K1vvdg4H1T8BnJ60O99/6YXebV6yfvAi4GnjgVr31/k84AROTjRpN/UEQuN+QHAVep6tVanT5+PXBBK+sFwGtM+I3At4hs/fPMo3yq6rtV9Qbz76XAuVvm0SKmTaE6lPdLwI3bZK6FGF6fDrxUzelxXfCQ4poQw6dSfUAJ4JZUX9PbOlT1vcBnBpJcAPyOVrgUOFtE7rgd7poY41VVL9H61oApx1RMu0J1lf6bgCn66Ch2wQV0vqo+QOuTc+cAn/TirzE0Hy6Nqs6pPlZzm00z2seDQYhPHz9EZWVNgVFezbT/Tqr6tm0yFkBMu94LuJeI/KWIXCoij94adzVi+Hwe8FQRuYbKAnz2dlhbGIv25V3BlGNqFCJyDtU5p9+cmpc+xBwES1gRIvJU4Dzg4VPzEoKIZMCvUt3PtA+YUbmBHkFlAb5XRO6rqp+blKsuvhd4taq+SEQeArxWRO6jqum7GitCRM6nUgAPm5qXAbwY+ClVLbfvoIjDpNtAReRjwGeppsq/paoXmYHyPFV9lEnzB1TT7X/Kyb/+LNLxg4SEhIRFcD2f/bQGLoObegbwMFW9VkRuD7xDRD4CXALc03yK8lrgHsCjVPXKW8it9RvlW6bkNyEhYRexoxb2ruCd5e9/IkSfVAGo6rXm9zoReTPwIFV9r4g8C3g71Qr6q1T1ShF5/s251ZTsJiQkJEF7pDCZAhCRmwKZql5vwt8GPB9Aq89DNj4Rqao/fwu59X/ZPqcJRwpJgCUkOEw5A/gK4M1mcWQGvE5V/3Q0VxrACQkJCWvBZApAVa8GJjvFl5CQMBFkF3afJ8D0i8CLI3WehISEhLVgvxSAgGTJBZSQkJCwEIoweb8UQKUBpmYiISEh4UhgVAGIyFnAc4A7q+rTReSewFep6jTf8U0zgISEhIS1IGYG8NvA3wAPMf9fC/w+E3zIXYBdPVKdkJCQsG+IUQD3UNXvEZHvBVDVGya4ebOCCOT5JFVvBBNew5GQkJAQowDOiMhNMN8HFpF7AKc3ytUQsrQGkJCQkLAOxCiA5wJ/CtxJRH4X+L+Y6tZIESRPCiAhISFhHRhVAKr6DhF5P9WXggT4MVX99MY5C0GA2Z5tXFo3ktsoISFhTeiVpoHvgn7K/N5ZRO6squ/fHFt9OGJrAAkJCQkTYsicfpH5PUX1MZO/pbLB7wdcTr0raHsQQQ4OVisjWdAJCQkJwIACUNXzwX2Q5YGq+nfm//tQfepu+xCBWZoBrAVJESbsI9I28LUixqH+VVb4A6jqh0TkazbIUz8E9OCYrwEsiiMi6KU8Gs8xJfSoHKJMSmBtiJGmV4jIK4D/af5/CnDF5lgagAicWNEFtA+YSmhPrSwGhLymid9mMZVymFqYT13/xIhRAD8APBP4MfP/e5noK/eaCeWpHZsBbOHz3rJuwbyu8tbN1w58Kn3tbR0BnVIIrXtX9TqfZc3tspV23rNd6jHbQG8Efs38TYtMKE8tMAOYzJJePmtIAEUXt4CbRBblcZm2XLb9V1QEqwjxfXM0rSzUlhVYi9S7II+6zCMtOYNZuf2m0t1rUmYxl8F9jMC4UNW7r4WDBaAiFKc24AtYcdQvLXAis42WPyYwI/gbVQhDZYzkHee/Fd/zimOVlku2SzOddVuf6xaqkQJ0UGCOFTHCc5Tg7ysjUpFFC/wFX9faZhdbVigx/pTzvPAp4EnArTfDzggymN9k/QpgKQEem2VUOA7l7Y8M5uuZAQTTBmjBdohNN6CENJC+wVPuEg4UEm7KwXcXOZMYVSxbdNFEW78DArtXGA2VHcgT5CVUdkD4BnkYqD82fZCnnrYYbMuBdzr6DlaN7613+9OJGBfQv7ZILxaRvwF+fjMsDfCSCfOz1uhkoVkvPQAAFK1JREFUc8IjruHHreQF844KV+lPFxJuXrpGORF1N/jz8nb4LtXx1czTLm/kOfz4sp0uUL8QVBAaeHeubisYFlGMnQo24BQaGOiWp1hBFy1ovXSNstv5veHVKLuVbKzeRh0tvvvq7zxz8Dni+APQkKgIKpUVlM9A+ui8wfK2owxiXED+ieCMakYwyUqsZnB41qJzs6HIkPBYvrx4C3qoXOnGNwS7T28mbMZJmG7K6/Dq11F2BbZLrxJWHh1exAneoKII6AdXR+kJBQJlDD2HJ+xdnoYyDCken+/QO1yjEgha2oF+6AsvG1+qC9vXK2hXwCLDAsyvLwvwIC1aSLBnEhDYIVqgvFaZYIR1SMl0aK2yQ3TDy1iesbhgub3pFneNLa4UFkwfgRhB/iIvPAc+Bnz3+lkZx1IKwCDazbOMy8YT3nHp+uvqs6olIDil3YN6BPtgXg3V6wld7aavlYGJK3G92cap+rS6XGfhll69Nt7SMh1239h6i7Ir5LXOKz7NKRcvzo+H5kyhHdcOL4qG8JPGr0AtxL04J1QyL18rLxm1wnRxitpLE4eGS1YLLndGQJpCHsz/LcUTFOziWd1+nAvXcYPKo52+TVuHohhSHq344P+x+frSjeQJYRMuohgF8EOqerVPEJG7rZ2TCGgG87OWzd0SRlEVDpQ2VE60uyccPzQDCFvlPi1kvbfy9ioFS/PayikAoxTKVnyrPHdgK0gT57pqpGvTSq1pvnAumuk0z6Aom7yUIKXJXFieS5euFvYlYmkufQlFVYk6WlHnMXGUWq9rlAFNZa4sF5FasNs7rLxvWoi92jzPwQpsQ9M8qz9+5NJlqC07NwKZzK2huFefZ7XAdPVLR9hrRi3kfVpLsGtOpzxfEAdpAWEfVDIDwrkh7MeUQofWM9tYdkYRyNtLGyojIl90GWtAjAJ4I9C+GO6NwNevn51haAbFqW3X6jMQnzTOv7xCWQpBpRZ098TROrOBhqIYmAGEZgqlrwC88lqCvVIozXTiKwCf1lIAUnphJ+y1ps09hWIVhYsrnZC3eZkXTqA3lENbeajWSiY0K/AFd9tiz7OGkHfpzBUnmtd5dWbSWQGbZ7XAntWCu473abTytoS8+W3TfCu+Vhh4wt7GdWkhwd6wzgOCPTSjsAhb+11ar+BexhofyLvJspbGimUP3Qb61cDXArcUke/yom5BtRtoZYjIo4Ffp7JhXqGqvziYIYPiJmv0x66ATWvmTWCp9Y0WLbTOMObjb8w2Ai4lly7kFvKUg2/l23S+cqnSaYDWCg/lbdcbUowQFvxthPzQQSHZtI7BCOeQgG27WLLhvHgCu523aZ37fHd5HrK2w1Z83ZCLLu6GsI/jzWLhMzdbxNAM4KuAbwfOBr7Do18PPH3VikUkB14KPBK4BrhMRP5IVT/cl0cFylUvAx3rSENvqyfvWnyIfr1DeRuDTQfjuzyNpI8pT7S7jtnDu4TaMsDfpjc89MrqAF1jJM0iAzqiuNh2qtIuUPcSCLbVWDsFjIJGeYG1qsG8gzPa0EAKh0NGxuA6XKDswdl1Dy12x99QvaPlrQlDt4H+IfCHIvIQVf2rDdT9IOAqu74gIq8HLgB6FQCA5v2tsdKULHa62H4boxZRQDgOletPcQdo1XpgU2CLqBeu02cuXZ03y7TOQ/W/SJOWi5JllZnsXMlZ2YgHmGUluUk3M2a1iDJr0WZZ0aEdZEUjXMUVHBh/j01/4NHcbzZ34dyMsAOZe+nmAJyQgszUccLEZZSNcPU8pSvHps+paRZZcA9uGGVrL2KBUBiz3MYVCIUJlybujOaNsE1vw4c6M7+5K+/QxeUclrMuzYTnZZ1+buNLE6dZI2zTzwM0qwQsrSgzV7ZtsaLMKEr7vKZNyozS5LW/qkJZ1uEqndQKyUunLmzjPIXkpW/n9WeggzQPDYXYnr026uvJ01NuU4w0Xaq0omLk/7JKYsgF9J9V9ZeB77MfhG/wpfqjy1XpcA7wSe//a4BvDPBxIXAhQH6rWy0/FRy0jLukXus8lGdI8I9Y08Ey2nk9mi/YG4Ifqq2DNkwdJy1hXymAslFelpVOoDsBn6kT7L7Qt0LZF/ou7GhdYT/LCg68MMBJT4j7wv5UdujC9vdki3ZC5pySdrq5E+xWARxIwQlaysMT9panAxTngjftlwO5aSRflOcRpn3hDV2rMgpVu5bt0eDQlHfoKQUbtoL7DLkn0Kuhe0bzhjIAuFEPONOinS4PGsoA4EaP5iuF00Z5WKF/mIUVgFUeuclbSEae1coATNc1fccK6UKgMMI+szQVt+htFYGIrxRajWYLB+MaMHXURMSE1b0H/50N0agFgHhRTmeYcaIyLFdcFd1yVQJiRehK+x72OlUtKReHXEB/b34vX67o9UBVLwIuAjh5pzupFAs+aaxqlP5/Rw+KBC36doEhN0mgA0lr1mDTB2aztkwNpOtVLiYu6MZpKRQk4HIQbSiXdnnSVlQNWjddI96rJmspo3a4SkMHIXdKFtsH9gBl0Ert0tpn3/w02rK+wTNkPQu7To9ndUuD7tN8d4+jQUc6qWd1Nyz3IfdRwNIOWvHtuAatm05URtxC/WWOunpi3UMD6Zvxy1q//RhyAf2x+X3N2mutcC1wJ+//cw2tF6JgjMCeBNHExbBEEd13OVLIMmyuvz9UGO2I/cm0E6jDfYvFQ9tPQ4uxgwu0SndxN0BrblP14lye7kOO+oPbaMwIu/TGQStLC+yyCS7u+unFy2PiOtssR2naW2+77NDzuXSt/yWQblPDNMjABvIsxeo6bZA1ljXkAvrjoapU9fEr1n0ZcE9zpuBa4MnA9w3mKCE/vUaptyHDMNrgjBS0o2UPlLPIwbKQMA2mjxK6WocDu3sa2zwdTd1vWzg3tnwObA0lsDU0K+qtoYS2i9ptoGXZCLu49qGwwvNDDG0DhXp/v3/Ay9vqCVT/+2FobO/Eo5U27NHcTp/BLZ/S2P4JNHYfhbeBGlrfjqT2SeFQ3qDia4WH0tGNjzKEY9OPlRUpbra+S2lN9Q25gF64nirCUNW5iDwLeDuVq/VVqnrlUB5RmH15Uwytnm7R6V2vRdmZpobjhvby12k8QbyMYB+yuv2rHlp5K8HeEtjqCWqbrvCsbV+wu/iuwF5KiPuHvWy6edGhuYNgljafezSbvqjjzZFmLRVx2zHtYa2sPvRlD4LleUUHZDZzNP+wF4DOujSybCXl0TiwZdK7dN7ef/fxHSf0JaAAxLkqB88BZD0nfEMKoCXUxhRF6KqKhZXH2IwmJGhj043kWSrNmjHkAnqPDYvICeCrqYbtR1X1zDoqV9WLgYtj00sJsxvWUXObkeXjl7r/Z9H4Ecs9nlYL35g8lWDXTlw3Xdfar9K1BHuovNCpXz9voZ28tfIIXAXRWCSsLWKhZYlrhlqhHLgyolZUdQM1L9iLsBhC1z/4fPnXP3SueOjSNBc6izIlrjHFrS7Xyki9BVUnpKxSmIsnML2ZQEtRBK9pyLy8A7MC9daWhsprxnv/d5SCNMsmkKdFa5Q/pADafPTxN5p/QJKvacYRizHFFHMZ3OOAlwH/QMXe3UTkGar6J+tgcBFUCiDWVI9AZFErWfaD6XQkvj9/KK5J027aBRSEo4VcO22eNPDM3r0/Q5fCNWcoXoKhnZbOgs26gj8Tj6+6PA0oiO69SYH0ISzywZqhy2v7FAQtIRK6GK59Nw94gluGBYl/55J9TFusUt/T5MrwDnO5OjTAM94UMPAsTmGEBLcEaARoft46waBCISKuj+ee+OFy+vvOSlvVN4DYy+DOV9WrAETkHsDbgEkUwME6FYCPtSiD/sh4pRCbLqw8xq6JHso7rGQCCsXCP0Eb4mlA0fUK/SHfutZCxu1JMlatqHpViOPPhe2VPBpo1ugLAxfog9EfIOmmCwqLmCuhG4IzJLWsMqzj7FZhvKbyy7Mkpxyg0xFUuvxVOyXt+woIbFdHVylU6fqfraFIYvK24vt5gdCgCwv78fe7FvdQBJZZh4hRANdb4W9wNdVp4O2jVGZf3syHYxc421MhWmEMJOyLGsgSVhBdYjBd4F78pRRObLqgMorlYbF2G2znTDrPHhosEjsSN3AcN2rw9nz4xN3A6hRko2QvnQ20FCm4d9UnsF0ZAaFr6xW/vjHXSKw1HVR4gXTBd2LbJZQ+QIuJo61YhpVQDILfLdgSYhTA5SJyMfB7VE/7JKprG74LQFX/YIP8NSAl61UAC3xDtw8LbTFfh9IYKSsqb18TRuTtfd5lFJ0rM4bn4TRLbfXfxMdelkRHOIdQhPmtzyyNPI8TpkMWRo9ws+SCDlyOHt7bZ6r662thRDCu4mvf3KchF0vfwZLfNl4WMQrgFPDPwMPN//8C3ITqfiAFtqgAlPyG+baqC/OwDpmxjOBZMMvCn7ncxAflF+VhFYG8mYnhBFi8DSTWggwI77YAHhU/C85+xg7LOiwg+JozjggsKlNXnOFtfUvoCoj5JOQPbIORKJSQnw714gUwpdW37rrXIPSW/qC9j3U/17be0S4qjWXcAasMiXW6s0IuowWyx1rlg6nW7U7Z9O17E9cbswvobsCzgbv66ddwEGxhiCpyeugo8IawQ4JiLQI7FttWljvkkpkMK9o3S2Pbgm5FZbEIpvjY+igm9Pv7iHEBvQV4JfDHTC0KVZHTE7iAjopg2qfn2Cde9xm7KByHsAS/e/aEFXZlBgDcqKov2TgnMSgVOTPBDOCoIQnXhKOGfVNkO4IYBfDrIvJc4M+A05aoqu/fGFd9UIXDaReBEyZCUlq7gSRojxRiFMB9ge8HvpnaBaTm/y1Dq7tZEhIYOa2bsBOQpDB2GjEK4EnA3dd1/89KUEWTAkhI2BskFb3biFEAH6L6LvB1G+ZlHErzKt6EhGWhqR8ltBB9oOLoIEYBnA18REQuo14DUFW9YHNs9UHd1bwJCX1IrqGE5XD8jIIYBfBcLyzAv6H6eMv2odR3sickJCQkrISYk8DvEZGvo/pa15OAj1FdD719aJoBJCQkJKwLQ5+EvBfwvebv08AbAFHV87fEWwdKmt4nJCQkrAtDM4CPAP8L+HbvWwA/vhWuhrCGGzwTEhISEoYVwHdR+frfLSJ/CryeyU9Va9q9kZCQkLAmDH0T+C3AW0TkpsAFwH8Ebi8ivwm8WVX/bEs8ekyBphlAQkJCwloQswj8JeB1wOtE5FZUC8E/RXU1xPaRZgAJ+4hjuMc8YfcRsw3UQVU/C1xk/hISEmKRDJeEHcRCCmAnkHYBHV+ke2USEtaKSRSAiDwPeDrV5yUBflZVL56Cl4Q9QlL+qyEp0IQWppwB/JqqvnDC+hMSjhf2UYEmpbVR7J8LKCEh4fhgU0orKRagOtm7/UorF9DTgC8AlwPPMQvMobQXAheaf+9DdTvpPuC2VCeodx37wifsD6/7wickXjeBXeTzLqp6uzZxYwpARN4J3CEQ9XPApVQNpMALgDuq6g9GlHm5qp63VkY3hH3hdV/4hP3hdV/4hMTrJrAvfMIGXUCq+q0x6UTk5cBbN8VHQkJCQkIYk5xOEZE7ev9+J/vj1klISEg4MphqEfiXReQBVC6gjwPPiMy3TwfQ9oXXfeET9ofXfeETEq+bwL7wOc0icEJCQkLC9EgXlCQkJCQcUyQFkJCQkHBMsZMKQEQeLSIfFZGrROSnA/EnReQNJv59InLX7XMZxedPiMiHReQKEflzEbnLFHwaXgZ59dI9QURURCbbxhbDq4h8t2nbK0Xkddvm0fAw9v7vLCLvFpEPmD7w2In4fJWIXCciwc0WUuEl5jmuEJEHbptHj5cxXp9iePw7EblERO6/bR49XgZ59dJ9g4jMReSJ2+ItGqq6U39ADvwDcHfgBPC3wL1baX4EeJkJPxl4w47yeT5wlgk/cwo+Y3k16W4OvJfqnMZ5u8orcE/gA8CtzP+331E+LwKeacL3Bj4+UZt+E/BA4EM98Y8F/oTqg08PBt43BZ+RvD7Ue++P2WVevX7yLuBi4IlT8dr3t4szgAcBV6nq1ap6hupLZBe00lwAvMaE3wh8i8jWz3aP8qmq71bVG8y/lwLnbplHi5g2hepQ3i8BN26TuRZieH068FI1p8dV9bot8whxfCpwCxO+JfCPW+SvZkL1vcBnBpJcAPyOVrgUOLu1VXtrGONVVS/R+taAKcdUTLsCPBt4EzBFHx3FLiqAc4BPev9fY2jBNKo6Bz4P3GYr3AV4MAjx6eOHqKysKTDKq5n230lV37ZNxgKIadd7AfcSkb8UkUtF5NFb465GDJ/PA54qItdQWYDP3g5rC2PRvrwrmHJMjUJEzqE65/SbU/PSh3QZ3BYgIk8FzgMePjUvIYhIBvwq1f1M+4AZlRvoEVQW4HtF5L6q+rlJuerie4FXq+qLROQhwGtF5D6q6eswq0JEzqdSAA+bmpcBvBj4KVUtt++giMMuKoBrgTt5/59raKE014jIjGp6/a/bYa/Dg0WIT0TkW6nuP3q4qp7eEm9tjPF6c6qL9v7CdNQ7AH8kIo9X1cu3xmWFmHa9hsr3ewh8TET+XyqFcNl2WATi+Pwh4NEAqvpXInKK6qKwXXMHRPXlXYGI3A94BfAYVd32uF8E5wGvN2PqtsBjRWSu1ffWdwNTL0IEFk1mwNXA3agX1762leY/0FwE/r0d5fPrqBYK77nrbdpK/xdMtwgc066PBl5jwrelcl/cZgf5/BPgaSb8NVRrADJRu96V/oXVx9FcBP7rKXiM5PXOwFXAQ6fkMYbXVrpXs4OLwDs3A1DVuYg8C3g71Qr6q1T1ShF5PnC5qv4R8Eqq6fRVVIswT95RPn8FuBnw+8YK+D+q+vgd5XUnEMnr24FvE5EPAwXwk7plSzCSz+cALxeRH6daEH6aGmmwTYjI/0PlLrutWY94LnBgnuNlVOsTj6USrDcAP7BtHi0ieP15qvW+/2HG1FwnunkzgtedR7oKIiEhIeGYYhd3ASUkJCQkbAFJASQkJCQcUyQFkJCQkHBMkRRAQkJCwjFFUgAJCQkJxxRJASQkJCQcUyQFkJCQkHBM8f8DRCQJ+XCvt2MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make prediction\n",
    "prediction = model.predict(X_star, operator=None)\n",
    "\n",
    "u = griddata(X_star, prediction[:, 0], (X, T), method=\"cubic\")\n",
    "v = griddata(X_star, prediction[:, 1], (X, T), method=\"cubic\")\n",
    "\n",
    "h = np.sqrt(u ** 2 + v ** 2)\n",
    "\n",
    "\n",
    "# Plot predictions\n",
    "fig, ax = plt.subplots(3)\n",
    "\n",
    "ax[0].set_title(\"Results\")\n",
    "ax[0].set_ylabel(\"Real part\")\n",
    "ax[0].imshow(\n",
    "    u.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "ax[1].set_ylabel(\"Imaginary part\")\n",
    "ax[1].imshow(\n",
    "    v.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "ax[2].set_ylabel(\"Amplitude\")\n",
    "ax[2].imshow(\n",
    "    h.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Gxlw5jxzXD-L"
   },
   "source": [
    "**Assess the accuracy**\n",
    "\n",
    "We can employ the test data to assess the accuracy of the trained model.  \n",
    "Using the normalized $L^2$ norm defined as  \n",
    "  \n",
    "Error $u$ = $\\frac{1}{||u_{test}||_{L^2}} ||u_{test} - u_{pred}||_{L^2}$  \n",
    "  \n",
    "where $u_{test}$ is the reference solution and $u_{pred}$ is the prediction given by the model, we get:  \n",
    "  \n",
    "Error $u$: 1.854433e-03  \n",
    "Error $v$: 2.413796e-03  \n",
    "Error $h$: 1.426797e-03  \n",
    "  \n",
    "We can also plot the absolute value of the error for each point of the domain. It's everywhere in the order E-3.  \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "j2SxMTchnP6C"
   },
   "source": [
    "![Error.png]()"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "Schrodinger.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
